### =========================================================================
### Preparation and settings ####
### =========================================================================

### Please set the working directory

wkpath <- ""
setwd(wkpath)

### Please set the folder with the occurrence data. 
### Data should be in .csv format (or any format delimited by ","), 
### with two columns "x" and "y" to give the longitude and latitude of occurrences
input_dir <- "example/cleaning_cc_occurrences"

# Please load polygon mapping functions from where you save these scripts
source("./functions/range_fun.R")

### Please set the names for output maps for each species, and the length should be the same with input occurrence .csv files
### If this is not set, the names would be the same with the input files.
  
  species_name <- NA # (optional)

### Default parameter settings:
  
  proj <- CRS("+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs +towgs84=0,0,0")   # Spatial projection in which the coordinates of the occurrence data (input) are stored. The output raster of the species range will have the same projection.
  Climatic_layer <- NA   #(Optional) climate raster (e.g. temperature) used to improve the distribution range (by rejecting cells with unsuitable climate).
  degrees_outlier <- 5   # distance threshold (degrees) for outlier classification. If the nearest minimal distance to the next point is larger than this threshold, it will be considered as an outlier. 
  clustered_points_outlier <- 2   # maximum number of points which are closer to each other than the degrees_outlier, but should still be considered as outliers.
  final_resolution <- 0.2   # determines the final resolution of the species range raster.
  buffer_width_point <- 0.5   # buffer (in degrees) which will be applied around single observations.
  buffer_increment_point_line <- 0.5   # how much should the buffer be increased for each point on a line.
  buffer_width_polygon <- 0.1   # buffer (in degrees) which will be applied around distribution polygons (for each bioregion)
  cut_off <- 0.05   # quantile of temperatures (on each site) which are not considered to be suitable for the species. e.g: cut_off=0.05: lowest and highest 5% of temperature distribution
  method <- "PCM"   # PCM (percentage cells method) method uses temperature filter only at the edge of range raster, ACM (all cells method) method for all cells.
  cover_threshold <- 0.3   # only if method=F. Specifies the threshold of proportion of covered area by the range polygon above which the corresponding raster will be classified as "present" without considering the temperature.
  dest_output='3.1_mapping_Polygons'   # path to where all rasters, plots, log files should be saved (if write_plot/write_raster=T)


### Please run following codes to get the polygon maps in the folder "Polygons" under your working directory

### =========================================================================
### Preparations of packages, functions and data ####
### =========================================================================

### install packages globally

# set a local mirror server
options(repos=structure(c(CRAN="https://stat.ethz.ch/CRAN/")))

# Packages to load
packages <- list("rgbif", "geometry", "maps", "ClusterR", "mclust", "raster","rgdal" ,"maptools" ,"tools" ,"ecospat" ,"rgeos" ,"PresenceAbsence")

# Install missing ones and load all packages
for (p in packages) {
  if(!library(package = p, logical.return = TRUE, character.only = TRUE)){
    install.packages(p)
    library(package = p, character.only = TRUE)
  } else {   
    library(package = p, character.only = TRUE) 
  }
}


### Download data from TNC. The shapefile will be saved in the folder of this script.
### If failed, please download at https://geospatial.tnc.org/datasets/b1636d640ede4d6ca8f5e369f2dc368b/about
### and then unzip and load manually
if (!dir.exists("./terr-ecoregions-TNC/")) {
  dir.create("./terr-ecoregions-TNC/")
  
  download.file(url = "https://www.arcgis.com/sharing/rest/content/items/b1636d640ede4d6ca8f5e369f2dc368b/data",
                destfile = "terr-ecoregions-TNC.zip", mode = "wb")
  
  unzip(zipfile = "terr-ecoregions-TNC.zip", exdir = "./terr-ecoregions-TNC")
  
  cat(paste0("The TNC bioregion map is downloaded and saved in the folder: \n",  getwd(), "/terr-ecoregions-TNC"))
}

if (list.files("./terr-ecoregions-TNC/", pattern = ".shp")){
  print("Loading TNC bioregion map.")
  
  bioregion <- readOGR(dsn="./terr-ecoregions-TNC", layer="tnc_terr_ecoregions")
} else {
  cat("The TNC Bioregion map does not exist. \n Please load your own map, or download and unzip from: \n https://geospatial.tnc.org/datasets/b1636d640ede4d6ca8f5e369f2dc368b/about")
}

### =========================================================================
### Polygon mapping ####
### =========================================================================

### Load data

species_csv_all <- list.files(input_dir, pattern = "\\.csv", full.names = T, recursive = T)
if (is.na(species_name)){
species_name <- gsub("\\.csv","",list.files(input_dir, pattern = "\\.csv", recursive = T))}


### Load occurrence and map the range

for (indexi in 1:length(species_csv_all)){
  
  species_name_i <- species_name[indexi]
  species <- read.csv(species_csv_all[indexi])
  
  occ_coord <- unique(species[,c("x", "y")])
  
  print(paste0("working with        ", species_name_i, "         [index:", indexi, "]"))
  
  species_range <- range.fun(species_name_i, occ_coord, proj=proj,  
                             Climatic_layer=Climatic_layer, Bioreg=bioregion, Bioreg_name = 'WWF_REALM2', 
                             final_resolution=final_resolution, degrees_outlier=degrees_outlier, clustered_points_outlier=clustered_points_outlier,
                             buffer_width_point=buffer_width_point, buffer_increment_point_line=buffer_increment_point_line, buffer_width_polygon=buffer_width_polygon, cut_off=cut_off, 
                             method=method, cover_threshold=cover_threshold, desaggregate_points=F, 
                             write_plot=T, write_raster=T,  
                             return_raster=T, overwrite=T,dest_output=dest_output,
                             dir_temp="./temp")
  
  # Bioreg: shapefile containg different bioregions (convex hulls will be classified on a bioreg basis)
  # Bioreg_name: how is the slot containing the bioregion names called?
  # desaggreagte_points: should close points be desaggregated? Speeds up clustering
  
  tmp <- tempfile()
  do.call(file.remove, list(list.files("tmp", full.names = TRUE)))
}

cat(paste0("Polygon maps are saved in the folder: \n",  getwd(), "/",dest_output))